$title	GTAP8DATA.GMS	Read a GTAP 8 dataset

$if not set nd $set nd 0

scalar nd	Number of decimals /%nd%/;
abort$(nd<>round(nd)) "Number of decimals must be an integer";

$if not set ds $set ds gsd
$if not set yr $set yr 07
$if not set datadir $set datadir "..\data\"
$setglobal datadir %datadir%

$if not defined g   set	g(*)	Goods plus C and G;

$if exist %ds%.gdx $gdxin '%ds%.gdx'
$if not exist %ds%.gdx $gdxin '%datadir%%ds%.gdx'
$if not defined g $load g

set		f(*)	Factors;

$load f

$if not defined i set i(g)	Goods;
$if not defined i	$load i

	

$if declared hs6 $load hs6
$if not defined r	set r(*)	Regions;
$if not defined r	$load r


set	rnum(r)	Numeraire region,
	sf(f)	Sluggish primary factors (sector-specific),
	mf(f)	Mobile primary factors;

alias (r,s), (i,j);

display r;


parameters
	vfm(f,g,r)	Endowments - Firms' purchases at market prices,
	vdfm(i,g,r)	Intermediates - firms' domestic purchases at market prices,
	vifm(i,g,r)	Intermediates - firms' imports at market prices,
	vxmd(i,r,s)	Trade - bilateral exports at market prices,
	vst(i,r)	Trade - exports for international transportation
$if declared hs6	viws_hs6(i,hs6,r,s)	Bilateral trade at world prices (HS6 goods)
	vtwr(i,j,r,s)	Trade - Margins for international transportation at world prices;

$load vfm vdfm vifm vxmd vst vtwr 
$if declared hs6 $load viws_hs6 

if (nd>0,
	vfm(f,g,r) = vfm(f,g,r)$round(vfm(f,g,r),nd);
	vdfm(i,g,r) = vdfm(i,g,r)$round(vdfm(i,g,r),nd);
	vifm(i,g,r) = vifm(i,g,r)$round(vifm(i,g,r),nd);
	vxmd(i,r,s) = vxmd(i,r,s)$round(vxmd(i,r,s),nd);
	vst(i,r) = vst(i,r)$round(vst(i,r),nd);
	vtwr(i,j,r,s) = vtwr(i,j,r,s)$round(vtwr(i,j,r,s),nd);
$if declared hs6	viws_hs6(i,hs6,r,s) = viws_hs6(i,hs6,r,s)$round(viws_hs6(i,hs6,r,s),nd);
);


set tartype		tariff type (ad-valorem or specific) /adv, spe/;

parameter
	rto(*,r)	Output (or income) subsidy rates
	rtf(f,g,r)	Primary factor and commodity rates taxes 
	rtfd(i,g,r)	Firms domestic tax rates
	rtfi(i,g,r)	Firms' import tax rates
	rtxs(i,r,s)	Export subsidy rates
$if declared hs6	rtms_hs6(hs6,r,s)	Bilateral tariff rates
	rtms(i,r,s)	Import taxes rates
	rtmssa(tartype,i,r,s)	Import taxes rates;


$load rto rtf rtfd rtfi rtxs rtms rtmssa
$if declared hs6 $load rtms_hs6

if (nd>0,
	rto(g,r) = rto(g,r)$round(rto(g,r),nd);
	rtf(f,g,r) = rtf(f,g,r)$round(rtf(f,g,r),nd);
	rtfd(i,g,r) = rtfd(i,g,r)$round(rtfd(i,g,r),nd);
	rtfi(i,g,r) = rtfi(i,g,r)$round(rtfi(i,g,r),nd);
	rtxs(i,r,s) = rtxs(i,r,s)$round(rtxs(i,r,s),nd);
	rtms(i,r,s) = rtms(i,r,s)$round(rtms(i,r,s),nd);
	rtmssa(tartype,i,r,s) = rtmssa(tartype,i,r,s)$round(rtmssa(tartype,i,r,s),nd);

$if declared hs6  rtms_hs6(hs6,r,s) = rtms_hs6(hs6,r,s)$round(rtms_hs6(hs6,r,s),nd);
);
parameter
	esubd(i)	Elasticity of substitution (M versus D),
	esubva(g)	Elasticity of substitution between factors
	esubm(i)	Intra-import elasticity of substitution,
	etrae(f)	Elasticity of transformation,
	eta(i,r)	Income elasticity of demand,
	epsilon(i,r)	Own-price elasticity of demand;

$load esubd esubva esubm etrae eta epsilon

parameter
	evdd(i,g,r)		Volume of domestic energy demand (mtoe),
	evdi(i,g,r)		Volume of imported energy demand (mtoe),
	evt(i,r,r)		Volume of energy trade (mtoe)
	eco2d(i,g,r)		Emissions of CO2 from domestic source (Gg),
	eco2i(i,g,r)		Emissions of CO2 from imported source (Gg)
	eco2(i,g,r);

$load evdd evdi evt eco2d eco2i
if (nd>0,
	evdd(i,g,r) = evdd(i,g,r)$round(evdd(i,g,r),   nd);
	evdi(i,g,r) = evdi(i,g,r)$round(evdi(i,g,r),   nd);
	eco2d(i,g,r) = eco2d(i,g,r)$round(eco2d(i,g,r),   nd);
	eco2i(i,g,r) = eco2i(i,g,r)$round(eco2i(i,g,r),   nd);

	evt(i,r,s) = evt(i,r,s)$round(evt(i,r,s),   nd);
);

eco2(i,g,r) = eco2i(i,g,r) + eco2d(i,g,r);


*	Declare some intermediate arrays which are required to 
*	evaluate tax rates:

parameter	vdm(*,r)	Aggregate demand for domestic output,
		vom(*,r)	Total supply at market prices;

vdm(i,r) = sum(g, vdfm(i,g,r));
vom(i,r) = vdm(i,r) + sum(s, vxmd(i,r,s)) + vst(i,r);

parameter
	rtf0(f,g,r)	Primary factor and commodity rates taxes 
	rtfd0(i,g,r)	Firms domestic tax rates
	rtfi0(i,g,r)	Firms' import tax rates
	rtxs0(i,r,s)	Export subsidy rates
	rtms0(i,r,s)	Import taxes rates;

rtf0(f,g,r) = rtf(f,g,r);
rtfd0(i,g,r) = rtfd(i,g,r);
rtfi0(i,g,r) = rtfi(i,g,r);
rtxs0(i,r,s) = rtxs(i,r,s);
rtms0(i,r,s) = rtms(i,r,s);

parameter	pvxmd(i,s,r)	Import price (power of benchmark tariff)
		pvtwr(i,s,r)	Import price for transport services;

pvxmd(i,s,r) = (1+rtms0(i,s,r)) * (1-rtxs0(i,s,r));
pvtwr(i,s,r) = 1+rtms0(i,s,r);

parameter	
	vtw(j)		Aggregate international transportation services,
	vpm(r)		Aggregate private demand,
	vgm(r)		Aggregate public demand,
	vim(i,r)	Aggregate imports,
	evom(f,r)	Aggregate factor endowment at market prices,
	vb(*)		Current account balance;

vtw(j) = sum(r, vst(j,r));
vom("c",r) = sum(i, vdfm(i,"c",r)*(1+rtfd0(i,"c",r)) + vifm(i,"c",r)*(1+rtfi0(i,"c",r)));
vom("g",r) = sum(i, vdfm(i,"g",r)*(1+rtfd0(i,"g",r)) + vifm(i,"g",r)*(1+rtfi0(i,"g",r)));
vom("i",r) = sum(i, vdfm(i,"i",r)*(1+rtfd0(i,"i",r)) + vifm(i,"i",r)*(1+rtfi0(i,"i",r)));

vdm("c",r) = vom("c",r);
vdm("g",r) = vom("g",r);
vim(i,r) =  sum(g, vifm(i,g,r));
evom(f,r) = sum(g, vfm(f,g,r));
vb(r) = vom("c",r) + vom("g",r) + vom("i",r) 
	- sum(f, evom(f,r))
	- sum(j,  vom(j,r)*rto(j,r))
	- sum(g,  sum(i, vdfm(i,g,r)*rtfd(i,g,r) + vifm(i,g,r)*rtfi(i,g,r)))
	- sum(g,  sum(f, vfm(f,g,r)*rtf(f,g,r)))
	- sum((i,s), rtms(i,s,r) *  (vxmd(i,s,r) * (1-rtxs(i,s,r)) + sum(j,vtwr(j,i,s,r))))
	+ sum((i,s), rtxs(i,r,s) * vxmd(i,r,s));

vb("chksum") = sum(r, vb(r));
display vb;

*	Determine which factors are sector-specific 

mf(f) = yes$(1/etrae(f)=0);
sf(f) = yes$(1/etrae(f)>0);
display mf,sf;

parameter       mprofit Zero profit for m,
                yprofit Zero profit for y;

mprofit(i,r) = vim(i,r) - sum(s, pvxmd(i,s,r)*vxmd(i,s,r)+sum(j, vtwr(j,i,s,r))*pvtwr(i,s,r));
mprofit(i,r) = round(mprofit(i,r),5);
display mprofit;

yprofit(g,r) = vom(g,r)*(1-rto(g,r))-sum(i, vdfm(i,g,r)*(1+rtfd0(i,g,r))
        + vifm(i,g,r)*(1+rtfi0(i,g,r))) - sum(f, vfm(f,g,r)*(1+rtf0(f,g,r)));
yprofit(i,r) = round(yprofit(i,r),6)
display yprofit;

*	Define a numeraire region for denominating international
*	transfers:

rnum(r) = yes$(vom("c",r)=smax(s,vom("c",s)));
display rnum;


* -------------------------------------------------------------------
* extra parameters for PCIintrade project
* -------------------------------------------------------------------

* -- SOME DATA FILTERING :
* drop negligible trade flows
vxmd(i,r,s) = round( vxmd(i,r,s), 8);
rtms(i,r,s)$(vxmd(i,r,s) = 0) = 0;
*vtwr(i,r,s)$(vxmd(i,r,s) = 0) = 0;
*(i,r,s)$(vxmd(i,r,s) = 0) = 0;




* -- data to extract from gtap :

parameter
*       caplab          capital labor ratio
                fd              final demand
                domfd		domestic final demand
		impfd           imported final demand
		absorbtion	absorbtion
*                expend          final demand = expenditure = C
                expshare(*,*)   sectoral shares of total expenditure (sigma)
                labor(r)        value-added generated by labor
                income(r)       income
                gdp(r)          gdp
                expenditure(r)  expenditure
                pcGDP(*)        per capita gdp (m)
                pcEXP           per capita expenditure
                btrade          bilateral trade flows (CIF)
                btariff         bilateral tariff
                btcost          bilateral transport cost
		Endow(f,r)      endowment in factor f

;

alias(i,ii);

gdp(r)=sum(g,vom(g,r))-sum((i,g),vdfm(i,g,r)+vifm(i,g,r));

* 2007 POPULATION FIGURES ARE TAKEN FROM GTAP8 DIRECTLY
* this is taken from GSDDAT.GDX, POP PARAMETER

parameter pop(*) /

omn	2.72630095481873
aus	21.0766525268555
nzl	4.22830009460449
xoc	9.43824195861816
chn	1317.88500976563
hkg	6.92589998245239
jpn	127.770751953125
kor	48.4560012817383
mng	2.61145305633545
twn	22.8784503936768
xea	24.2407531738281
khm	14.323842048645
idn	224.66960144043
lao	6.0923318862915
mys	26.5556545257568
phl	88.7181854248047
sgp	4.58860015869141
tha	66.979362487793
vnm	85.154899597168
xse	50.5784378051758
bgd	157.752517700195
ind	1124.78698730469
npl	28.2867279052734
pak	162.590896606445
lka	20.0100002288818
xsa	29.2367324829102
can	32.976001739502
usa	301.290008544922
mex	105.280517578125
xna	0.126463994383812
arg	39.4904632568359
bol	9.52449512481689
bra	190.119995117188
chl	16.6361351013184
col	44.3594436645508
ecu	13.3418169021606
pry	6.12664318084717
per	28.5084800720215
ury	3.32390594482422
ven	27.4829998016357
xsm	1.50363600254059
cri	4.45878219604492
gtm	13.3537693023682
hnd	7.17412900924683
nic	5.59505176544189
pan	3.34334111213684
slv	6.10676097869873
xca	0.311500012874603
xcb	40.3928604125977
aut	8.30078792572021
bel	10.6256999969482
cyp	0.853814005851746
cze	10.3341598510742
dnk	5.46143817901611
est	1.34167194366455
fin	5.31641530990601
fra	63.5008583068848
deu	82.2663726806641
grc	11.1927633285522
hun	10.0557804107666
irl	4.35693120956421
ita	59.3752899169922
lva	2.27609992027283
ltu	3.37561798095703
lux	0.479992985725403
mlt	0.409049987792969
nld	16.3816967010498
pol	38.1205596923828
prt	10.6083345413208
svk	5.39731788635254
svn	2.01812195777893
esp	44.8789443969727
swe	9.14809226989746
gbr	61.0051116943359
che	7.55111694335938
nor	4.71127033233643
xef	0.346906989812851
alb	3.13245797157288
bgr	7.6597638130188
blr	9.70199966430664
hrv	4.43599987030029
rou	21.5468730926514
rus	142.100006103516
ukr	46.509349822998
xee	3.6674690246582
xer	14.2805843353271
kaz	15.4841995239258
kgz	5.23479986190796
xsu	38.5725631713867
arm	3.07244992256165
aze	8.58129978179932
geo	4.35785722732544
bhr	0.759559988975525
irn	71.0210418701172
isr	7.18009996414185
kwt	2.66296601295471
qat	1.13755297660828
sau	24.1574306488037
tur	73.0037384033203
are	4.36391305923462
xws	86.0129852294922
egy	80.0605392456055
mar	31.2241363525391
tun	10.2253999710083
xnf	40.3998413085938
cmr	18.6599388122559
civ	20.1227951049805
gha	22.8709659576416
nga	147.72184753418
sen	11.8933353424072
xwf	81.4210968017578
xcf	20.6531581878662
xac	80.0773696899414
eth	78.646125793457
ken	37.7546997070313
mdg	18.6043643951416
mwi	14.4394960403442
mus	1.2606920003891
moz	21.8693618774414
tza	41.2762107849121
uga	30.637544631958
zmb	12.3139419555664
zwe	12.44921875
xec	72.9722747802734
bwa	1.89242601394653
nam	2.0886709690094
zaf	47.850700378418
xsc	3.18307495117188
xtw	0.00656800018623471

/;

* transform to millions :
pop(r) = pop(r) * 1000000 ;




domfd(i,r) = vdfm(i,"c",r)+vdfm(i,"g",r)+vdfm(i,"i",r);
impfd(i,r) = vifm(i,"c",r)+vifm(i,"g",r)+vifm(i,"i",r);
fd(i,r)= domfd(i,r)+impfd(i,r);
expenditure(r) = sum(i, fd(i,r));

* test sensitivity to inclusion of g and i
parameter fd_disagg, fd_disagg_share;
fd_disagg("c",i,r) = vdfm(i,"c",r)+vifm(i,"c",r);
fd_disagg("g",i,r) = vdfm(i,"g",r)+vifm(i,"g",r);
fd_disagg("i",i,r) = vdfm(i,"i",r)+vifm(i,"i",r);

fd_disagg_share("c",i,r) = ( vdfm(i,"c",r)+vifm(i,"c",r)) / sum(i.local, vdfm(i,"c",r)+vifm(i,"c",r));
fd_disagg_share("g",i,r) = ( vdfm(i,"g",r)+vifm(i,"g",r)) / sum(i.local, vdfm(i,"g",r)+vifm(i,"g",r));
fd_disagg_share("i",i,r) = (vdfm(i,"i",r)+vifm(i,"i",r)) / sum(i.local, vdfm(i,"i",r)+vifm(i,"i",r));

execute_unload 'fd_disagg.gdx', fd_disagg, fd_disagg_share;



impfd(i,r) = vifm(i,"c",r)+vifm(i,"g",r)+vifm(i,"i",r);
fd(i,r)= domfd(i,r)+impfd(i,r);



* calculate per-capita values:
pcEXP(r)$pop(r) = 10e8*expenditure(r)/pop(r);
Endow(f,r) = sum(g,vfm(f,g,r));

btariff(i,s,r) = (1+rtms0(i,s,r)) * (1-rtxs0(i,s,r));
btrade(i,r,s) = btariff(i,r,s)*vxmd(i,r,s) + sum(j, vtwr(j,i,r,s))*pvtwr(i,r,s);
absorbtion(i,r) = sum(g, vdfm(i,g,r)+ vifm(i,g,r));


parameter taxrevenue(r); 

taxrevenue(r) =  sum(j,  vom(j,r)*rto(j,r))
	+ sum(g,  sum(i, vdfm(i,g,r)*rtfd(i,g,r) + vifm(i,g,r)*rtfi(i,g,r)))
	+ sum(g,  sum(f, vfm(f,g,r)*rtf(f,g,r)))
	+ sum((i,s), rtms(i,s,r) *  (vxmd(i,s,r) * (1-rtxs(i,s,r)) + sum(j,vtwr(j,i,s,r))))
	- sum((i,s), rtxs(i,r,s) * vxmd(i,r,s));


* -- LOAD DISTANCE DATA
* updated with set of GTAP 8 regions

set dist /       contig,comlang_off, colony,  distw, ldist,  homebias /;


parameter  importdist(r,s,dist);


*$if not exist importdist.gdx $call gdxxrw i=cepii_data.xlsx o=importdist_gtap8.gdx par=importdist rng="sheet1!a1:F50177" cdim=1 rdim=2
$gdxin %datadir%importdist_gtap8.gdx
$load importdist


*display importdist;

*fix Romania
* done this in excel file instead
*importdist("rom",r,dist) = importdist("rou",r,dist); 
*importdist(r,"rom",dist) = importdist(r,"rou",dist); 


importdist(r,s,"homebias")$sameas(r,s) = 1; 

set nodistdata;
nodistdata(r)$(not sum(s,importdist(r,s,"distw"))) = yes;
display nodistdata;
* ok, only the residual regions.

parameter regstats;

regstats("total") = card(r);
regstats("nodistdata") = card(nodistdata);

display regstats;

importdist(r,s,"ldist")$importdist(r,s,"distw") = log(importdist(r,s,"distw")); 

